The impact of agricultural diversification on U.S. crop production: National Models Estimation and Diagnostics
Kate Nelson
11/30/2018
This file contains models examining the functional relationship between agricultural landscape diversity and agricultural yields. This work applies hierachical Bayesian spatiotemporal modeling techniques to estimation of the effect of diversity of the agricultural landscape on the yield of corn, soy, and winter wheat.
REPRESENTING SPACE
The yield data available is at a county scale and the distribution of yields across space exhibits strong autocorrelation where yields in neighboring counties are more alike than yields in distant counties. This spatial autocorrelation is accounted for using a standard Conditional Autoreggressive dependency model based on adjacency for all counties in the conterminous US. In order to account for additional county-specific factors that contribute to yields a county iid random effect term is also included, yielding a Besag-York-Mollie (BYM) spatial dependency model. A county adjacency matrix is created for each crop under investigation. For regional models a seperate county adjacency matrix for each region-crop combination is created.
MODELS
The model form is a random-effects panel model with county spatial effects, non-parametric SDI, non-parametric Temp and Precip controls, and a quadratic time trend (similar to Schlenker and Roberts, 2009) for each agro-eco-region. The spatial effects account for standard, unstructured or iid, county random effects (deviation of the county mean from the population mean) as well as conditional autogregressive (CAR) structured residuals between counties. These effects account for time-invariant factors that differ across counties (however assume that the effects are not correlated with errors hence do not fully account for potential ommitted variables). The quadratic time trends for each agro-eco-region account for factors that are space-invariant (within an agro-eco-region) and vary across time (without the overfitting that was observed with the use of the rw2 random time effects). This effect is intended to account for factors that produce changes over time, but that may differ across agro-eco-regions. These effects are expected to partially account for differences in technology use and adoption rates over time in different regions and differences in agricultural management styles and techniques across regions that result in temporally varying responses to climate conditions. Note that while there clearly are differences in management styles, climate responses, and technology use within agro-eco-regions these effects account for aspects of the agro-eco-region as a whole that constrain or enable certain types of actions.
formula<-YIELD ~ 1 + f(TP, model=“rw1”,scale.model=T) + f(SDD,model=“rw1”,scale.model=T) + f(GDD,model=“rw1”, scale.model=T) + f(SDI, model=“rw1”, scale.model=T)+ PERC_IRR + ACRES + f(CNTY, model=“bym”, graph=H.adj.corn, scale.model=TRUE) + f(AERCODE.id, YEAR.id)+ f(AERCODE.id2, YEAR.id2)
Previous models employed r-inla default prior settings for the random effects and used a reduced precision fixed effect prior. This round of models keeps the reduced precision fixed effect prior and further employs penalized complexity priors (https://arxiv.org/pdf/1403.4630v4.pdf) that use a scaling factor to specify priors based on the limits of the data creating a weakly informative prior. Default and recomended settings for PC priors as provided in the R-INLA documentation and in https://arxiv.org/pdf/1403.4630v4.pdf are employed. These models also employ sum-to-zero constraints for all random effects(https://link.springer.com/article/10.1007%2Fs00477-017-1405-0).
ABOUT DIAGNOSTICS
In the returned model ouputs several diagnostics are reported. Some of these diagnostics are computed by R-INLA, and additional posterior predicve and cross-validation checks were manually computed. Descriptions of the diagnostics are adapted from Spatial and Spatio-temporal Bayesian Models with R-INLA (2015) by Blangiardo and Cameletti.
The deviance information criterion (DIC) is a generalization of the AIC (Akaike information criterion) developed for Bayesian model comaprison that examines tradeoffs between model fit (log-likelihood) and model complexity (effective number of parameters). Models with smaller DIC are better supported by the data.
The conditional predictive ordinate (CPO) is a leave-one-out-cross-validation diagnostic that estimates the probability of observing a value in the predictive distribution given a model fitted to all the observed data except that value. Larger values of CPO indicate a better model fit.
The probability integral transform (PIT) is another leave-one-out-cross-validation check that evaluates the predictive performance of the model where a Uniform distribution of PIT means the predictive distribution matches well with the data.
Manually calculated posterior predictive checks (PPC) were obtained using the predictive distribution estimated by R-INLA. These checks use all observations for both model evaluation and predictive checks (hence are not a form of cross-validation).
The first PPC shows a scatterplot of the posterior means for the predictive distribution and observed values, where points scattered along a line of slope 1 indicates a very good fit.
The second PPC provides a histogram of the posterior predictive p-value which estimates the probability of predicting a value more extreme than the observed value, where values near 0 and 1 indicate that the model does not adequately account for the tails of the observed distribution. The tails are due to isolated county-year events which seem to correspond with disaster records.
The final two PPC metrics calucalted provide the mean squared error (MSE) and R squared (R2) between the observed and predicted data (Blangiardo and Cameletti, 2015, page 168).
In addition, an independent cross-validation model was run using 6/7ths of the data for model fit estimation and reserving a held-out sample of 1/7th of the data for comparison. Three diagnostic metrics (MSE, R2, and NSE) were calculated for the cross-validation models. As a smaller sample size is used for these caluculations they are more sensitive to outliers than the previously reported diagnostics. Therefore, outlier error values greater than twice the mean of the observed data were removed prior to calculation of these diagnostics and the count of outliers is also reported.
The mean squared error (MSE) and R squared (R2) between the held-out data and predicted data were calculated (Blangiardo and Cameletti, 2015, page 168).
The Nash Sutcliffe Efficiency (NSE) which is commonly used to assess the predictive power of hydrological models was computed. A value of 0 indicates that the mean of the observed data is as good a predictor as the model, values less than 0 indicate the model is a poorer preditor than the mean of the data, and a value of 1 indicates that the model perfectly predicts the observed data.
US EXTENTS CORN MODELS
This section estimates the model described above for corn for three metrics of agricultural landscape diversity (SIDI, SDI, and RICH). Model estimates, diagnostics, functional response plots, and maps of spatial effects are all reported.
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ---------- ---------- ----------- ---------- -----------
(Intercept) 4.0998 0.0218 4.0569 4.0998 4.1425
PERC_IRR 0.0091 0.0006 0.0079 0.0091 0.0102
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 29.9498 0.4522 29.0791 29.9435 30.8537
Precision for TP 323.2940 153.4361 110.6326 295.8486 700.5702
Precision for SDD 20.8065 5.5533 11.9756 20.1122 33.6745
Precision for GDD 114.8603 46.6638 50.5940 105.9077 230.4463
Precision for SIDI 2003.0531 1482.3413 440.7825 1609.4662 5900.7110
Precision for CNTY 22.2986 1.1455 20.1984 22.2404 24.7028
Phi for CNTY 0.9990 0.0017 0.9941 0.9997 1.0000
Precision for AERCODE.id 355.2637 81.5375 221.4464 346.3688 540.6241
Precision for AERCODE.id2 9911.1632 2139.6328 6373.4412 9685.9031 14747.1260
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
--------- --------- ---------- ----------
-4756.95 2209.102 0.0289812 0.7082578[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.7857263 | 0.0382403 | 0.6607148 | 2 |
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ---------- ---------- ----------- ---------- -----------
(Intercept) 4.1323 0.0212 4.0908 4.1323 4.1740
PERC_IRR 0.0088 0.0006 0.0076 0.0088 0.0100
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 29.9944 0.4531 29.1064 29.9934 30.8910
Precision for TP 323.9264 157.6665 113.1573 292.9302 715.6133
Precision for SDD 20.6454 5.5220 11.8783 19.9501 33.4501
Precision for GDD 112.2088 45.1103 49.0859 103.9458 223.3186
Precision for SDI 885.7269 600.4967 230.4002 729.7607 2465.8757
Precision for CNTY 22.0985 1.1524 19.7891 22.1258 24.3032
Phi for CNTY 0.9993 0.0013 0.9955 0.9999 1.0000
Precision for AERCODE.id 371.9377 85.8291 231.0064 362.6058 567.0494
Precision for AERCODE.id2 9913.3866 2181.6194 6380.7361 9652.6657 14916.9165
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
---------- --------- ---------- ----------
-4772.191 2214.884 0.0289038 0.7084921[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.7397979 | 0.0423634 | 0.6614034 | 2 |
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ----------- ---------- ----------- ----------- -----------
(Intercept) 4.1761 0.0230 4.1310 4.1760 4.2213
PERC_IRR 0.0084 0.0006 0.0073 0.0084 0.0096
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 29.9331 0.4532 29.0635 29.9248 30.8437
Precision for TP 341.6731 174.5707 118.6400 304.0864 785.2311
Precision for SDD 21.8095 5.7672 12.5552 21.1180 35.1047
Precision for GDD 117.2635 47.1568 51.7564 108.4094 233.3232
Precision for RICH 1079.5893 566.2113 350.2677 960.4435 2510.8869
Precision for CNTY 22.4371 1.1727 20.0898 22.4638 24.6855
Phi for CNTY 0.9992 0.0014 0.9951 0.9999 1.0000
Precision for AERCODE.id 412.9991 97.1548 248.1269 404.8015 628.0932
Precision for AERCODE.id2 11123.1604 2410.7423 6922.2920 10965.9996 16340.8399
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
---------- --------- ---------- ----------
-4762.183 2213.329 0.0289657 0.7080809[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.7112438 | 0.0413223 | 0.6744548 | 2 |
US EXTENTS SOY MODELS
This section estimates the model described above for soy for three metrics of agricultural landscape diversity (SIDI, SDI, and RICH). Model estimates, diagnostics, functional response plots, and maps of spatial effects are all reported.
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ----------- ----------- ----------- ---------- -----------
(Intercept) 3.2781 0.0141 3.2504 3.2781 3.3057
PERC_IRR 0.0065 0.0005 0.0054 0.0065 0.0075
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 49.7800 0.7936 48.2088 49.7863 51.3274
Precision for TP 695.5598 493.8363 133.6313 574.7619 1977.3435
Precision for SDD 72.6201 21.8427 38.2643 69.8419 123.2019
Precision for GDD 287.6908 106.0468 123.2948 274.2293 533.0373
Precision for SIDI 15200.4452 22564.7069 1532.3233 8593.0572 69739.3076
Precision for CNTY 33.1739 1.7864 29.8960 33.0843 36.9240
Phi for CNTY 0.9974 0.0031 0.9889 0.9986 1.0000
Precision for AERCODE.id 250.6527 55.3843 158.1464 245.1963 375.1507
Precision for AERCODE.id2 7275.8801 1535.8317 4640.2622 7154.3799 10643.1035
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
---------- --------- --------- ---------
-9225.362 4393.392 0.017492 0.751808[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.7736001 | 0.0258817 | 0.7081603 | 2 |
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ---------- ---------- ----------- ---------- -----------
(Intercept) 3.2882 0.0131 3.2624 3.2882 3.3138
PERC_IRR 0.0064 0.0005 0.0054 0.0064 0.0074
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 49.8351 0.8359 48.1051 49.8774 51.3696
Precision for TP 663.5682 465.6018 123.1728 552.2777 1868.7598
Precision for SDD 72.3292 21.7891 37.3316 69.8613 122.2682
Precision for GDD 285.3265 105.9807 121.5044 271.7122 530.9557
Precision for SDI 2760.1632 1710.5156 626.3132 2390.6095 7070.1910
Precision for CNTY 33.5448 2.2422 28.8027 33.7247 37.4636
Phi for CNTY 0.9973 0.0036 0.9872 0.9988 1.0000
Precision for AERCODE.id 259.5931 61.5184 165.0557 250.4716 404.9463
Precision for AERCODE.id2 7220.0371 1633.9353 4664.8174 6992.7377 11034.4481
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
--------- --------- ---------- --------
-9241.69 4401.045 0.0174666 0.75212[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.733186 | 0.0242373 | 0.7366562 | 2 |
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ---------- ---------- ----------- ---------- -----------
(Intercept) 3.2933 0.0133 3.2672 3.2933 3.3194
PERC_IRR 0.0064 0.0005 0.0054 0.0064 0.0074
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 49.8778 0.7920 48.3326 49.8735 51.4507
Precision for TP 679.8154 524.4508 139.9951 538.2348 2061.3030
Precision for SDD 71.6777 21.6642 37.6102 68.9289 121.7865
Precision for GDD 284.4466 105.1302 125.8213 269.2925 533.5239
Precision for RICH 4902.6156 3740.9666 1252.4640 3843.1526 14820.7651
Precision for CNTY 33.0099 1.7383 29.7691 32.9423 36.6026
Phi for CNTY 0.9981 0.0024 0.9914 0.9991 1.0000
Precision for AERCODE.id 245.0916 56.0029 155.4249 238.0331 374.5956
Precision for AERCODE.id2 7383.8038 1557.0342 4647.7376 7289.4864 10732.7118
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
---------- --------- ---------- ----------
-9242.452 4399.599 0.0174389 0.7522154[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.7599772 | 0.0246851 | 0.7117944 | 2 |
US EXTENTS WINTER WHEAT MODELS
This section estimates the model described above for winter wheat for three metrics of agricultural landscape diversity (SIDI, SDI, and RICH). Model estimates, diagnostics, functional response plots, and maps of spatial effects are all reported.
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ---------- --------- ----------- ---------- -----------
(Intercept) 3.7795 0.0297 3.7217 3.7793 3.8384
PERC_IRR 0.0069 0.0006 0.0057 0.0069 0.0081
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 30.3264 0.5634 29.1978 30.3369 31.4122
Precision for TP 194.1723 82.3368 74.5400 181.1520 392.3227
Precision for SDD 659.1264 310.7045 257.3424 592.3727 1448.5368
Precision for GDD 531.0536 277.8829 178.1044 470.7707 1239.1956
Precision for SIDI 926.9797 648.8700 229.7683 756.4284 2633.5675
Precision for CNTY 12.4511 0.8780 10.7602 12.4427 14.2149
Phi for CNTY 0.9756 0.0109 0.9499 0.9772 0.9919
Precision for AERCODE.id 134.2762 26.6501 89.7140 131.6345 194.0044
Precision for AERCODE.id2 4486.9209 898.7138 3039.6699 4374.9881 6548.7653
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
---------- --------- ---------- ----------
-3410.903 1585.825 0.0281662 0.7687149[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.7669969 | 0.0404798 | 0.7480352 | 2 |
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ---------- --------- ----------- ---------- -----------
(Intercept) 3.8526 0.0317 3.7911 3.8523 3.9153
PERC_IRR 0.0063 0.0006 0.0051 0.0063 0.0075
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 30.2602 0.5574 29.1759 30.2563 31.3698
Precision for TP 186.7372 82.9091 75.7365 170.0672 393.9371
Precision for SDD 615.3965 278.2376 241.1038 560.4547 1310.6071
Precision for GDD 501.0939 254.4968 174.8109 446.5999 1147.2450
Precision for SDI 403.7058 324.3212 90.2708 312.1551 1261.6378
Precision for CNTY 12.5667 0.8843 10.9004 12.5447 14.3715
Phi for CNTY 0.9750 0.0110 0.9488 0.9766 0.9912
Precision for AERCODE.id 140.7010 27.9104 93.6204 138.0694 203.0532
Precision for AERCODE.id2 4584.8992 906.0762 3101.9769 4480.3011 6653.0673
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
---------- --------- ---------- ----------
-3397.041 1578.491 0.0281994 0.7682754[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.7362974 | 0.0391073 | 0.7371566 | 2 |
[[1]]
Table: Summary Table of Model Estimates
mean sd 0.025quant 0.5quant 0.975quant
---------------------------------------- ---------- --------- ----------- ---------- -----------
(Intercept) 3.8501 0.0300 3.7917 3.8499 3.9095
PERC_IRR 0.0063 0.0006 0.0051 0.0063 0.0075
ACRES 0.0000 0.0000 0.0000 0.0000 0.0000
Precision for the Gaussian observations 30.4295 0.5616 29.3509 30.4205 31.5553
Precision for TP 198.8608 88.1498 79.5291 181.6452 418.3981
Precision for SDD 660.4567 308.0338 256.7287 596.0708 1436.6173
Precision for GDD 530.9166 290.7141 181.0598 461.9682 1283.2551
Precision for RICH 568.7480 257.2415 218.6604 519.3383 1208.6633
Precision for CNTY 12.6909 0.9147 10.9541 12.6714 14.5458
Phi for CNTY 0.9731 0.0116 0.9459 0.9748 0.9904
Precision for AERCODE.id 136.0891 26.7152 90.2826 133.8813 194.9322
Precision for AERCODE.id2 4616.9970 879.1737 3107.1747 4544.2102 6555.3250
[[2]]
[[3]]
Table: Model Diagnostic Metrics
DIC CPO MSE R2
---------- --------- ---------- ----------
-3437.564 1598.147 0.0280755 0.7692436[[1]]
[[2]]
[[1]]
[[1]]
[[1]]
[[1]]
| R2 | MSE | NSE | OutlierCount |
|---|---|---|---|
| 0.7940621 | 0.0404092 | 0.7583557 | 2 |
SUMMARY PLOTS
Summary plots showing the effect of the smooth, random-walk, effects of diversity and climate on log(Yield) for each model. All crops are presented on a single plot, and new plots are created for different diversity metrics and different climate variables. Solid lines represent the median effect and the colored bands represent the 95% credibility limits as obtained from the estimated posterior distribution.
Effect of Diversity Metrics
Effect of Climate Variables
